Shear Induced Structural Ordering of a Model Metallic Glass 
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We report results of non-equilibrium molecular dynamics simulations of a one-component glassy 
system under the influence of a shear flow, with the aim of investigating shear induced ordering of 
this system. In spite of the very low temperature, the system transforms into a strained crystalline 
state through well defined nucleation events. Various characteristics of the observed ordering at 
different strain rates and temperatures are discussed. We also define and discuss the transition 
rates. 
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Nucleation, clustering and ordering processes in super- 
cooled melts are a subject of considerable interest, both 
in soft matter and in hard materials such as metallic al- 
loys or molecular glass formers. Crystal ordering is gener- 
ally thought to proceed through nucleation events, which 
are described within the framework of classical nucleation 
theory (^CNT)^^^^. In this theory, the nucleation is as- 
sociated with the local crossing of a free energy barrier 
under the influence of thermal noise. The situation be- 
comes more intriguing if, instead of a supercooled liquid, 
we consider a glassy system at very deep supercooling^. 
The kinetic slowing down and high viscosity of a glass 
make the structural rearrangements, including those re- 
quired for crystal nucleation, so slow that the correspond- 
ing transition is practically unobservable on the experi- 
mental time scales. Obviously, the devitrification can be 
accelerated through reheating of the glass. An interest- 
ing alternative, which is in fact commonly used in soft 
matter systems, is to utilize shear as an ordering tool 
rather than temperature. Crystal-crystal or amorphous 
crystal transformations are also known to take place in 
processes such as grinding and milling of molecular ma- 
terials, which in general involve a complex combination 
of mechanical stresses and temperature variations^. 

Recent numerical studies provide the evidence of 
crystallization induced by external shear in an amor- 
phous, single-component Lennard- Jones system^^^ How- 
ever, such single-component Lennard- Jones systems have 
a very weak glass forming ability and the transition to a 
crystalline state at low temperature proceeds through a 
spinodal, rather than nucleation, route^^^^. 

In this study, we consider the case of a system with 
a much better glass forming ability than the Lennard- 
Jones system, but still keeping the simplicity associated 
with one-component systems. This system is made of 
particles interacting through the short ranged, oscillating 
Dzugutov potentialiiii^ 



/7(r*)/e = A (r*"^-5) exp 



— a 



e(a-r*) 



B exp 



6(6 -r*), r* = r/a, (1) 



where a and e are the length and energy scales, respec- 
tively, the parameters A = 5.82, B = 1.28, m = 16, 
a = 1.87, b = 1.94, c = 1.1, d = 0.27 are taken as de- 
scribed in RefJ^ and 9(r) represents the Heaviside step 
function. So, the potential is characterized by a mini- 
mum U{rmin) = — 0.581e located at Vmin = 1.13cr and 
a maximum U{rmax) = 0.46e at Tmax — 1.628<j. The 
presence of these oscillations can be considered as a sim- 
plified attempt to reproduce the screened Coulomb re- 
pulsion of ions in metals. The potential is tuned in such 
a way that it prevents the formation of the most com- 
mon crystal phases at low pressure. At high pressures 
(P > 5e/cr^) and low enough temperatures^^, the FCC 
or HCP crystalline structures can be formed, although 
the ground state of the system at zero temperature is a 
BCC structure. A slow quench from a high-temperature 
state this model favors the formation of quasi-crystals in- 
cluding Frank-Kasper structures, whereas a faster cool- 
ing leads to a metastable amorphous state with a mode- 
coupling temperature Tmct estimated to be OAe/kB 
(see Ref.^^). A remarkable property is that this is one- 
component good glass former shares many common fea- 
tures with much more complex metallic glasses. With 
only one component, however, the study of topological 
reorganization, which excludes such effects as chemical 
ordering or demixing, is facilitated^^-. 

The system consists of = 19 652 particles con- 
tained in a cubic simulation box of the volume with 
L = 28.55(7. The structural changes in such a rather 
large system is expected to be relatively insensitive to the 
boundary conditions imposed by the simulation cell^^. A 
set of well equilibrated independent samples is created at 
a temperature T = 2.0e//cB. These samples are cooled at 
a rate AT/ At = 0.000 97e/ (kBr) and the configurations 
corresponding to the temperatures T = 0.01, 0.03 and 
0.06e//cB are stored^^. To impose the shear we create 
two solid walls at the ends of the simulation box in the 
y direction by freezing all the particles in the x-z plane 
over the range of three interparticle distances. The top 
wall moves in the x direction with the instantaneous ve- 
locity v{t) = jLy{t)ex at a constant strain rate 7 and 
normal load Pyy = 7.62e/cr^; Ly is the distance between 
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the walls. This value of the normal pressure is the 
namically favorable to FCC crystalline forms anc 
one to avoid the phase of quasi-crystals for the 
system^^. Constant temperature conditions are < 
by rescaling the velocity component in the neutr 
rection perpendicular to the flow x and the veloc 
dient y direction. The termostatting was imple 
only if the temperature has changed more than 
the considered value. Other details of the simulat: 
tocol are identical to those described in Refj^. 

The local order in the system can be charactei 
parameters that probe the appearance of orien' 
order around each atom. In this work, we use b 
local orientational^ bond-order parameter qi{i) 
global^ counterpart Qf. 
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Here Yi^{Oij ^(fij) are the spherical harmonics, Nij{i) de- 
notes the number of the nearest "neighbors" of particle 
i, and Oij and cpij are the polar and azimuthal angles 
formed by radius-vector r^j and some reference system. 
The parameter Qi is invariant under rotations. Nonzero 
values of Qi for the different crystalline cluster struc- 
tures including those with icosahedral symmetry appear 
for / > 6 (see Ref.^^). In this work we define "neighbors" 
as all atoms located within a sphere of radius |r^j | < 1.5cr 
around an atom i. For a fully disordered system one 
has Qe ^ 0, whereas the ordered geometries are char- 
acterized by the next values of Qq: Q^^^ = 0.5745, 
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0.4848 and 



^6 - 0.5106, QI""' = 0.6633, 
= 0.3536 (see RehM^). 
By performing a cluster analysis for the detection of 
the "solidlike" particles (for details see Refs.— ^2^), we ob- 
tain the cluster size distribution ns{t) and its moments 
s^^\t) = s'^ns{t) as a function of time. The ze- 
roth moment s^^^ defines the total number of clusters, 
the first moment 5^^^ represents the total number of par- 
ticles involved in an ordered phase (solidlike particles), 
the ratio of the second and first moments, s^^^s^^^ is 
associated with the average cluster size. In a finite sam- 
ple we also obtain s* = max^^ n^^o[^s] as the size of the 
largest cluster. The average spatial cluster size {() can be 
also estimated through the distance between maximally- 
separated particles within the same cluster. 

We begin the discussion of our results with Fig. [H 
which represents the time evolution of the orientational 
order parameter Qeit)^ of the number of solidlike parti- 
cles s^^\t) and of the number of clusters s^^\t). These 
curves are computed from a single run in a sample 
sheared with the rate 7 = 0.005r~^ at the tempera- 
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Figure 1: Time-dependence of different characteristics com- 
puted for a single sheared sample with the strain rate 7 — 
0.005r~^ at the temperature T = 0.03e//cs: (a) Global ori- 
entational bond-order parameter Qeit). Inset: radial distri- 
bution of particles formative the ordered cluster structures. 
It is defined for the moment indicated by arrow. Pronounced 
peaks of radial distribution appear at the distances ~ cm^^^, 
n = 1, 2, 3, . . . is the number of a peak, that is typical 
for FCC structures, (b) Number of solidlike particles s^'^\t). 
Inset shows the same in logarithmic scale, (c) Number of 
clusters s^^\t). 



ture T = 0.03e/kB' All these quantities are correlated 
and illustrate clearly the transition the initially amor- 
phous system into an ordered state, which takes place 
rapidly after a system-dependent incubation time. Fur- 
ther, the transition is also seen from the growth in the 
average cluster size presented in Fig. [2l The clear tran- 
sition from an amorphous phase to a strong crystalline 
order after an incubation time, during which the system 
stays amorphous, distinguishes this case from Lennard- 
Jonesium^. Such a behaviour, in which the incubation 
time is a statistical, sample-dependent quantity, is typ- 
ical for a nucleation process. After the transition, the 
system is transformed essentially into a FCC structure, 
as is clearly evident from the snapshots and from the ra- 
dial distribution of particles constitutive of the ordered 
structures, gci{r) (inset of Fig. [Ifa)). Both the cluster 
size distribution and the structure of the clusters are es- 
sentially independent of shear rate after the transition is 
completed, indicating that the local structure is not shear 
rate dependent. It should be noticed that ordering takes 
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nucleus as a function of time is shown in Figs. [4land[5l 
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Figure 2: Time evolution of two average cluster characteris- 
tics in the sheared system with the strain rate 7 = 0.005t~^ 
at the temperature T = 0.03e/A;s. Main: Average cluster 
size as the ratio between the second and first moments of the 
cluster size distribution. Inset: spatial correlation length (see 
text). Relatively insignificant values of both terms for the 
ordered phase is consequence of the retaining of small crys- 
talline grains (for the range t > 600r the ratio between the 
largest cluster s* presented in Fig. [5] and the average cluster 
size s^^-'/s^"'"'* is ~ 4 times). 
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Figure 3: (Color online) Density profiles of solidlike particles 
for the instantaneous configurations at the different times af- 
ter shear initiation. The results are typical for all studied 
samples. Profiles presented here are for the same system as 
in the previous figure. 



place primarily away from the walls, as is seen from the 
density profiles of solidlike particles presented in Fig. [3l 
Similar profiles were observed for all considered samples 
under shear and for all considered strain rates and tem- 
peratures. Such a behavior points to the homogeneous 
scenario of the nucleation. The evolution of the largest 
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Figure 5: (Color online) Main: Size of the largest cluster 
versus time. Broken line corresponds to the value 160. The 
study sample is still the same as for the previous figures. Inset: 
Three largest clusters at the moment indicated by arrow. The 
sizes of the clusters are A^i = 177, A^2 = 130 and N3 = 126 
particles. 

After the transition into an ordered phase, which is 
accompanied by the growth of amount of solidlike parti- 
cles s^^^ [see Fig. [TJb)], the total number of clusters s^^^ 
[see Fig. [U^c)] as well as the size of the largest cluster 
presented in Fig. [5] become essentially stationary. The 
system achieves a high level of structural ordering, which 
is, however, far from perfect. When ^ 83 % of the bulk 
becomes included into an ordered phase, no further or- 
dering is observed. Remaining particles formative the 
disordered ranges are located mainly near the walls. Al- 
though the statistical properties remain essentially con- 
stant, the size of the largest cluster displays considerable 
fluctuations (Fig. [5]). The system is made of nanocrys- 
talline grains containing at most a few hundred particles 
(see inset of Fig.[5j), which are prevented by the shear flow 
from growing further but still exchange particles under 
the influence of the mechanical noise generated by the 
shear flow. 

Until now we have discussed the strain history of a sin- 
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Figure 6: (Color online) (a) Orientational bond-order pa- 
rameter of a sheared system (amorphous initially) versus 
strain. Each curve corresponds to a fixed strain rate 7 G 
[0.0001, 0.01]r-^ and is averaged over results obtained for at 
least 50 independent samples, which were sheared under the 
same conditions at a temperature T = 0.03e/A;s. (b) Main: 
Orientational bond-order parameter versus rescaled transi- 
tion strains ^ = l/ln- Inset: Mean first passage time dis- 
tributions multiplied by the corresponding strain rates. The 
solid lines are results of simulations, the dotted lines present 
the fits to Eq. {|3}. 



gle, typical sample. We now turn to the consideration of 
the statistical properties of this strain history by collect- 
ing statistics on 50 and more independent runs. The evo- 
lution of orientational bond-order parameter Q^i^) (with 
7 = 7t) at the temperature T = 0.03e//cB and at various 
values of the strain rate 7 is presented in Fig.[6](a). These 
curves display a marked transition at a well defined value 
7n of the strain, which increases with shear rate. Obvi- 
ously, the transition can be characterized by a transition 
rate r~^, which we evaluate using the mean first-passage 
time (MFPT) method. This method^^^^^^^^ focuses on 
the average time r needed for the system to reach for the 
first time some state with a defined value of the reaction 



coordinate ^ ^ 0i, ^2, ^3, G M, where R is the 
domain of states^^. As a rule, this allows one to evaluate 
the distribution of time r{(j)) over the reaction coordinate 
(j). For an activated process, the transition rate is 
simply determined by the inverse of the MFPT, at which 
a plateau is reached, i.e. = r{(j)p)~^^ whereas the 

total r{(j)) can be fitted by 

r{cj^) = ^{l+evf{{cj>-r)c)), (3) 

where erf(. . .) is the ordinary error function, c is the lo- 
cal curvature around the top of the barrier and 0* cor- 
responds to the critical value of reaction coordinate (for 
details, see Ref."*^^). 

The transition rates calculated within the framework 
of the MFPT approach with the order parameter as the 
reaction coordinate, Qq (/), were used to rescale the 
strain dependencies presented in Fig. [6Ka). So, Fig. [6l^b) 
shows the order parameter [same with Fig. [D^a)] with 
the rescaled argument ^ = 7/7n = t/rn^ where 7^ = jrn 
is the transition strain. As can be seen, all curves col- 
lapse onto a single master curve, which is indicative of 
a common mechanism for the transition scenario. The 
transition takes place in three stages, which, in turn, re- 
flect the influence of the external drive on the inherent 
dynamics of the system. The first initial stage, where the 
system remains mainly in an amorphous phase, is char- 
acterized by shear-induced "unjamming" of the system 
and by the fluidization of the glass^^. Due to the shear 
the diffusivity grows at this stage and the potential en- 
ergy increases. The characteristic time scale associated 
with inherent dynamics is defined here by the relaxation 
time at a finite shear rate, ra{j) (see Ref.^). The second 
stage corresponds directly to the transition of the amor- 
phous system into an ordered state and is caused by clus- 
ter nucleation^. Its appearance can be characterized by 
the transition time scale (or the transition strain 7^). 
During this stage the order parameter grows rapidly, fol- 
lowing a power law oc ^^^o ^. In the third stage, shear 
initiates rearrangements in the ordered structures and 
prevents the system from undergoing a complete crystal- 
lization, the order parameter increases only very slowly 
and probably saturates for large strains2ii2S. 

The extracted values of the transition rates at dif- 
ferent strain rates and three temperatures are presented 
in Fig. [71 At small values of 7, the transition rates in- 
crease with strain rate in a manner close to a power law. 
Note that the rates are practically independent of tem- 
perature at these low values of 7. However, with the 
subsequent increase of 7 the growth of transition rates 
becomes slower and temperature dependent. Eventually, 
starting with 7 ^ 0.005r~^ transition rates decrease with 
the strain rate 7. 

To describe the transition of a system into an ordered 
phase, the ordinary practice is to use the Arrhenius-like 
equation 

r-i=^oexp(-/3A^), (4) 



5 



fiuence of the shear. As a result, Eq. (|4j) can be rewritten 
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Figure 7: (Color online) Transition rate versus strain rate at 
three different temperatures. The points with error bars are 
the data extracted from our molecular dynamics simulations. 
The solid curves are the fits using Eqs. dH and (|6]). The short 
straight line corresponds to a unit slope. 



where is the attempt frequency (or kinetic prefactor) 
to overcome a barrier and (5 = l/{kBT) is the in- 
verse thermal energy; G = PAE. Two time scales arise 
naturally in this equation, and Uq^. Identifying a rate 
with the nucleation rate^^ and the barrier AE with 
the free energy cost to form a cluster of the critical size, 
Eq. dH) is the main equation of CNT. The use of Eq. 
(|4]) for systems near thermal equilibrium is straightfor- 
ward: from the measured transition rates at different 
temperatures one can easily define the thermodynamic 
barrier AE as well as the kinetic term uq. However, the 
use of this equation for non-equilibrium cases that vio- 
late detailed balance is not generally justified, although 
there are some indications that a similar approach using 
an "effective temperature" could be relevan t^^i^^ . 

Nucleation events in systems under shear have been 
identified in Ref*^ with the conclusion that the thermo- 
dynamic barrier is an increasing function of shear 
rate. We propose an interpretation of our results along 
the same lines by making a number of reasonable assump- 
tions, namely: 

(i) A formula similar to (|3j) can be used providing an 
unknown effective temperature Teff, which is introduced 
to replace the usual bath temperature T = l/{kB/3). 

(ii) The kinetic prefactor is set by the time scale for relax- 
ation in a system under shear, which is inversely propor- 
tional to the shear rate in this low temperature regime^. 
Hence uq cx |7|. 

While the assumption (i) results in the generalization 
of the energy barrier height G = [3AE in Eq. (|4]) to the 
7-dependent ternt^i^, G{j^T) = AE / (kBT^f^)^ assump- 
tion (ii) is motivated by the fact that the kinetic processes 
that can be observed in a glass are dominated by the in- 



a|7|exp[-G(7,T)]. 



(5) 



This is justified by the fact that the crystal nucleation in a 
glassy system free from any drive, e.g. 7^0, becomes at 
very deep supercooling extremely rare event, i.e. 
0, practically unobservable within the time-scale of the 
experiment as structural rearrangements are inhibited by 
the very high viscosity^. 

These assumptions are not sufficient to determine inde- 
pendently the effective temperature and the barrier AE^ 
both of which, as can be expected, depend on shear rate 
and on the actual, thermodynamic temperature. How- 
ever, the variation of a thermodynamic factor G(7, T) 
with shear can be defined. A first clear observation from 
Fig. [71 is that G{j,T) must be typically of order unity. If 
we make the reasonable assumption that the barrier AE 
is of the same order of magnitude as in usual supercooled 
liquids (typically of order e, where e is the energy scale of 
the interaction), it appears that the usual activation at 
the temperatures under study (T < 0.06e//cB) is unable 
to account for our observations. Moreover, any effective 
temperature Teff, which would be used to explain the ob- 
served values of G(7,T), should depend on shear rate 
only weakly. 
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Figure 8: (Color online) Change of the exponent of thermo- 
dynamic factor G{j,T) with the strain rate in regard to the 
case with 7 = O.OOOlr"^ 

The shear rate dependence of G(7, T) is presented in 
Fig. [8l As the effective temperature can only increase 
with shear rate, the data clearly shows that the energy 
barrier must increase with shear rate. This is consistent 
with the findings of Ref.^. However, we found that the 
parabolic fit (in the variable 7) of G(7,T), which was 
used in Ref*^ to treat the data, is not sufficient to de- 
scribe our results. Instead we use a cubic dependence of 
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Table I: Numerical data for different temperatures T on the 
parameter A and the correction factors Bi , B2 and Bs to the 
equilibrium barrier height G^^^ as obtained from the fitted 
simulation data. 



T {e/ks) 


A 


Bi (r) 


B2 (r^) 




0.01 


2.44 


430 


-230 X 10^ 


500 X 10^ 


0.03 


2.78 


417 


-229 X 10^ 


505 X 10^ 


0.06 


2.86 


390 


-210 X 10^ 


510 X 10^ 



shear rate: 

3 

G(7,T) = G(°) + ^i?„(r)|7r, (6) 

n=l 

where G^^^ characterizes the limit of the energy barrier 
hight at zero shear rate and B^iT) are the correction 
factors. Equation (j6|) is an expansion in powers of the 
dimensionless strain rate 7r as a small parameter, 7r G 
[0.00001, 0.01], for G(7,T) about the equilibrium value 
G^^\ Then, Eq. ([5]) can be extended to the following 
form: 

r-i=^|7|exp(^-^B„(r)|7r^, (7) 

where A oc exp [— G^^^] characterizes the zero shear rate 
limit of the thermodynamic factor. The fits of Eq. ([7|) 
are shown as continuous lines in Fig. [3 and the obtained 
numerical values of all parameters involved in Eq. ([7j) are 
given in Table [H Note that Eq. ([7j) should be considered 
as an attempt to extend the main equation of CNT for 
non-equilibrium case^ of a glassy system under shear 
at very deep supercooling, where the inherent structural 
rearrangements are extremely slow and the external shear 
affects the transport properties directly as was discussed 
in Refs.^^^^^. The fact that the coefficients Bn(T) are, in 
fact, almost independent of the temperature shows that 
the process is activated by shear, rather than by thermal 
fluctuations. 

In order to obtain an independent determination of 
the effective temperature, it would be necessary to de- 
termine an energy barrier Eq at zero temperature and 



zero shear rate for the formation of a crystalline nucleus, 
and to obtain Teff as the ratio Eo/G^^\ Unfortunately 
no method is known to obtain Eq directly. A possibil- 
ity, which we have not pursued in this work, would be to 
extrapolate results from high temperatures. As a result, 
our expression for the nucleation rate involving an effec- 
tive temperature cannot be seen as a direct validation of 
this concept, but rather as a plausible extrapolation of 
other results^Si^a. 

Our results demonstrate that the formation of a 
nanocrystalline phase in a driven, low temperature sys- 
tem proceeds via a nucleation mechanism, not dissimilar 
to the one observed in the standard liquids at moder- 
ate undercooling. Nucleation theory in its classical form 
is not applicable, in view of the very low temperatures. 
However, our results are consistent with a nucleation ap- 
proach that would involve an effective temperature and 
a nucleation barrier, which increases with shear rate, as 
already reported in supercooled liquids. 

This work could be extended in several directions, in 
particular: 

(i) Metallic glasses often have quasi-crystalline phases. 
For the Dzugutov system studied here, the quasi- 
crystalline phase can coexist with the liquid (or amor- 
phous) and crystalline phases (see Fig. 4 in Ref.^^), and 
the spontaneous growth of quasi-crystals in this system 
was recently investigated by Glotzer^^. It could be in- 
teresting to study the competition between quasi-crystal 
formation and crystal nucleation under an external drive. 

(ii) Glassy materials are characterized by a waiting 
time dependence of their properties, or aging, and a state 
that depends on the cooling rate. The cooling process 
and the waiting time could have an influence on the or- 
dering of the system under an external driving. 

(iii) Although the one-component system considered 
in this study is a good glass-former, bulk metallic glasses 
are, as a rule, multicomponent (five different kinds of 
atoms, at least). Therefore, it would be of interest to 
study the influence of structural inhomogeneity appeared 
from polydispersity on the induced structural ordering. 
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